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Abstract 

Aims. In this paper we present two efficient implementations of the diffusion approximation to be employed in Monte Carlo com- 
putations of radiative transfer in dusty media of massive circumstellar disks. The aim is to improve the accuracy of the computed 
temperature structure and to decrease the computation time. The accuracy, efficiency and applicability of the methods in various 
corners of parameter space are investigated. The effects of using these methods on the vertical structure of the circumstellar disk as 
obtained from hydrostatic equilibrium computations are also addressed. 

Methods. Two methods are presented. First, an energy diffusion approximation is used to improve the accuracy of the temperature 
structure in highly obscured regions of the disk, where photon counts are low. Second, a modified random walk approximation is 
employed to decrease the computation time. This modified random walk ensures that the photons that end up in the high-density 
regions can quickly escape to the lower density regions, while the energy deposited by these photons in the disk is still computed 
accurately. A new radiative transfer code, MCMax, is presented in which both these diffusion approximations are implemented. These 
can be used simultaneously to increase both computational speed and decrease statistical noise. 

Results. We conclude that the diffusion approximations allow for fast and accurate computations of the temperature structure, vertical 
disk structure and observables of very optically thick circumstellar disks. 
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1. Introduction 

Protoplanetary disks are the sites of planet formation and thus 
are the ideal environments for studying how planetary systems 
are formed. With the current development of new observational 
techniques and the increasing sensitivity of modern telescopes, 
the call for accurate, fast and flexible methods of interpreting 
these observations is increasing. Analysis tools are needed that 
are fast and, at the same time, that do justice to the rich detail of 
current observations. 

An important step in translating predictions from disk prop- 
erties into comparison with observations lies in solving the ra- 
diative transfer in protoplanetary disks subject to conservation 
of energy and thermal equilibrium. The opacity and energy bal- 
ance in protoplanetary disks are dominated by their dust grains, 
even though they only contribute about one percent of the total 
disk mass. Unfortunately, determining the radiation transport is 
far from a trivial task, mainly because of the wide range of op- 
tical depths that have to be traced in the disk. On the one hand, 
the optically thin upper layers of the disk are important since 
they are the regions where solid-state resonances are formed that 
give information on the size distribution and composition of the 
dust grains in these layers. On the other hand, the cold midplane 
region of the disk is important since it contains most of the mass 
of the disk and thereby determines the vertical structure. Radial 
optical depths through this midplane can easily reach values of 
10 5 -10 6 in the visible part of the spectrum. 

The required flexibility for radiative transfer can be obtained 



by using a Monte Carlo method (see e.g.JBjorkman & Wood 



is that it is accurate at all optical depths and puts virtually no 
restrictions on the spatial distribution of matter or on the opti- 
cal properties of the dust grains. The disadvantage is that the 
method becomes increasingly slow for high optical depths since 
every interaction of a photon package has to be considered sepa- 
rately. Also, in the deep layers of the disk that are well-shielded 
from the stellar radiation, photon counts may be low, causing 
large errors on the derived temperature structures in these re- 
gions. Accurate values for the temperatures in the midplane re- 
gions of protoplanetary disks are crucial for computing self- 



2001 ; |Niccolini et all [2003) |Pinte et al.| |2006| >. Its advantage 
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consistent vertical density distributions (Dullemond & Dominik 
[2004l|Dullemond et al.||2007) . 

In this paper we discuss solutions to these problems in terms 
of two different formulations of the diffusion approximation for 
high optical depth regions. The two different formulations serve 
different purposes and can be used together. We present a par- 
tial diffusion approximation (PDA) which can be used to de- 
crease statistical noise on the temperature structure and a modi- 
fied random walk (MRW) procedure which is useful for decreas- 
ing computation time significantly. Both used together provide 
a fast and accurate radiative transfer scheme that can be used 
for a variety of applications. We focus on passive protoplanetary 
disks, i.e. where the radiation from the central star dominates, 
and consider only the radiative transfer through the dust com- 
ponent. In Sect. [2] we will outline the general radiative transfer 
method employed along with the implementations of the PDA 
and MRW. Then, in Sect. [3] we will check the accuracy of the 
approximations and discuss the resulting self-consistent vertical 
density structures. The results are put together into recommenda- 
tions for the use of the diffusion approximation in various cases 
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in Sect. |4] Finally, in Sect. [5] we summarize our findings and 
draw our conclusions. 



2. Computational approach 

2.1. Monte Carlo radiative transfer 

A flexible way of computing the transfer of radiation in an in- 
homogeneous, complex medium is by means of a Monte Carlo 
method. In this method photon packages emerging from the cen- 
tral star are traced through the disk, allowing them to undergo 
scattering, absorption and reemission events caused by the dust 
they encounter on their path. The Monte Carlo method for ra- 
diative transfer in dusty media has become particularly efficient 
since the ideas of continuous absorption and reemission have 
been proposed and worked out by Bjorkman & Wood ( 2001| l. 
The method is fully described in their paper, and we will not 
repeat it here. 

Two major difficulties of Monte Carlo radiative transfer are: 

1. The statistical noise in regions shielded from radiation. 

2. The fact that computation time increases with roughly the 
square of the optical depth. 

At first glance, the first difficulty might seem not so rele- 
vant since in radiative transfer in protoplanetary disks the noise 
mainly appears in regions which are optically well shielded, i.e. 
mainly the midplane of the disk, and the midplane contributes 
only little to the observed infrared spectrum. However, for im- 
ages, visibilities and millimeter observations it can be of impor- 
tance to have a reliable estimate of the midplane temperature. 
In addition, the midplane temperature is crucial for hydrostatic 
equilibrium computations of the vertical structure of the disk. 
For this reason, computations of the vertical structure of proto- 
planetary disks have, up to now, relied on more simplified meth- 
ods that do not have the problem of noise. We will in the next 
section outline a method to circumvent this problem, and obtain 
a reliable approximation of the temperature in regions with low 
photon counts. 

2.2. Diffusion approximation for optically thick regions 

2.2.1. Decreasing photon noise: Partial Diffusion 
Approximation (PDA) 

In regions in the disk where the mean free path of the photons 
is much smaller than the length scale over which density, p, and 
temperature, T, change, the energy transport is by radiative dif- 
fusion. This is described by the radiative diffusion equation (see 
|Wehrse et al.| ( [2000] > who provide an extension of the classical 
one dimensional result by |Rosseland| ( |1924| to three dimensions) 

1 dE 

V-(DV£) = - — , (1) 
c ot 

where E is the local energy density, c is the speed of light, and D 
is the diffusion coefficient which in the Rosseland approximation 
is given by 

D R = U R = -4- (2) 
3 ipxR 

1 is the Rosseland mean free path, and^ is 



where A R = (px R ) 
the Rosseland mean opacity given by 
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Hereby is the wavelength dependent mass extinction coefficient, 
B v is the Planck function, and v is the frequency of radiation. 
The Rosseland diffusion coefficient is an approximate measure. 
In case scattering cannot be neglected, the diffusion coefficient 
has to be adjusted. The expression for the diffusion coefficient 
we use is given in section '. 



2.2.3 



We can use the diffusion equation to improve the accuracy 
of the temperatures in regions of high optical depth. The way we 
implement this is that after the Monte Carlo procedure we solve 
the time independent version of Eq. (QJ, i.e. we put dE/dt = 
in regions of low photon counts. In this case we get the diffusion 
approximation for the temperature in the cell 



V • (DVr 4 ) = 0. 



(4) 



-dv 



Eq. Q results in a system of linear equations which can be 
solved if we set appropriate boundary conditions. Since we wish 
to solve the diffusion equation only in a limited region of the disk 
(the regions around the midplane), the boundary conditions are 
simply set by the temperatures as determined by the Monte Carlo 
procedure at the edge of the region with low photon counts. 
We refer to this method as the Partial Diffusion Approximation 
(PDA). The accuracy of the results of this approximation will be 
discussed in detail later in this paper. Note that the PDA is only 
used after the full Monte Carlo run has already finished. Thus 
observables that are directly obtained from the escaping Monte 
Carlo photon packages are not influenced by it. The main use 
of the PDA is to obtain a reliable temperature structure for, for 
example, iteratively solving the vertical hydrostatic density dis- 
tribution. 

The PDA assumes that the probability that a photon escapes 
from the optically thick region without interactions is zero. In re- 
ality though, a photon emitted at a very long wavelength encoun- 
ters a very low opacity, and thus the disk can become optically 
thin for these photons. This will cause these regions to cool more 
efficiently than predicted by the PDA. Thus, in general the PDA 
is expected to overestimate the temperature slightly. We will re- 
turn to this point when we discuss the accuracy of the PDA. 



2.2.2. Increasing computational speed: Modified Random 
Walk (MRW) 

Monte Carlo radiative transfer computations can be very slow 
when implemented in a straightforward manner, especially for 
geometrical configurations where a photon package has a large 
number of interactions with the dust before it leaves the sys- 
tem. For protoplanetary disks this occurs especially when a sin- 
gle photon package 'gets lost' in the optically thick regions of 
the disk, i.e. the disk midplane. It can then take a very high num- 
ber of interactions before the photon package finally leaves the 
system. Although only very few photons get to these regions, in 
some configurations of rather massive protoplanetary disks these 
photons completely dominate the computation time. We can try 
to overcome this problem by letting photons make multiple inter- 
action steps in a single computation. This can be done by using 
the radiative diffusion equation as outlined below. 

In the regions of radiative diffusion the photons travel by a 
random walk. This random walk is basically a solution to the ra- 
diative diffusion equation (Eq.[T]i. In the Monte Carlo procedure 
we wish to know the distance a photon travels before leaving a 
given region in the disk. For the Monte Carlo procedure, space 
is subdivided in regions of constant density and temperature and 
we can apply the method of Fleck & Canfield ([1984 ). The solu- 
tion to Eq. ([T} in the case of an infinite homogeneous medium is 
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simply 



1 



{AnDctfl 2 



exp 



ADct 



(5) 



where if/{r, t) is the fraction of the energy that has diffused to 
position r in a time f, and r — |r|. 

Eq. ([5]) is only valid for an infinite, homogeneous medium 
while the cells in the disk are of finite size. The way to solve 
this is outlined by [Fleck & Canfie"Id| ([T984) and consists of set- 
ting an absorbing boundary condition to the diffusion equation 
at r — Ro, with Rq being the radial boundary of the homogeneous 
region. This ensures that the time it takes for the photon to reach 
the boundary of the region is the time when the photon crosses 
the boundary for the first time. The solution to this is a Modified 
Random Walk (MRW) and is given by 

To determine the likelihood that the photon is still inside a sphere 
with radius Rq after a time t we have to integrate iff(r, t) from 
to R 



Pit) = 4tt r 2 t/r(r, t)dr = 2 Y (-1)' 
Jo r~i 



with 



y = exp < -D ct 
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(7) 



(8) 



In the Monte Carlo procedure we now draw a random number, 
Po, between and 1 and solve Pit) = Po for t. This now gives us 
the time the photon has traveled before reaching the edge of the 
sphere. The energy that is absorbed by the dust grains per unit 
mass during this random walk to the edge of the sphere is equal 
to 

E = E y ctpk, (9) 

where E y is the energy of the photon package, and k is the aver- 
age mass absorption coefficient encountered by the photons 
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where r\ v is the emission coefficient, and k v is the wavelength de- 
pendent mass absorption coefficient. The emission coefficient r\ v 
is proportional to the probability that a photon leaves an interac- 
tion event (scattering or re-emission) with frequency v. An iter- 



ative way to obtain r\ v is given in section 2.2.3 The above equa 



tion averages the frequency-dependent mass absorption coeffi- 
cient weighted with the average distance traveled at that particu- 
lar frequency. This procedure can then be repeated until the pho- 
ton is close enough to the cell's edge where the regular Monte 
Carlo method can take over. 

Since the MRW is only valid in optically very thick regions, 
an important aspect of the above method is a criterion to start 
the MRW procedure. A good criterion for this is that the small- 
est distance to the cell's edge is larger than a few times the 
Rosseland mean free path (Fl eck & Canfield[|1984[ ), i.e. 



y 

pxr 



(ii) 



The parameter y can be adjusted for speed or accuracy. Values of 
y close to unity ensure fast radiative transfer but with relatively 
large errors, while higher values of y ensure more accurate re- 
sults but at the cost of a slightly increased runtime. In addition 
to this criterion we also require the photon package to have ex- 
perienced several tens of interactions in a single cell before the 
MRW is employed. Note that there is always a finite chance that 
a photon gets emitted at a wavelength where the optical depth is 
low and the photon can escape the cell without interactions. The 
criterion set in Eq. (Hi has to be chosen such that the fraction of 
such photons is small. We study the effects of different values of 
y in section [3~2| 



2.2.3. The Diffusion Coefficient 

In this section we derive the diffusion coefficient used in the 
modified random walk procedure and the partial diffusion ap- 
proximation. We start with the diffusion coefficient as given by 
|Fleck & Canfield| ( fT984| ) for the modified random walk 



D 



6(d) 



(12) 



where d is the distance a photon package travels, (d) is the mean 
free path, and (d 2 ^ is the mean of the pathlength squared. For 
monochromatic radiative transport at frequency v we have that 
(d) = d = ipXvT 1 and (d 2 ) = d 2 = 2ip Xv )- 2 . 

Since we have absorption and reemission in addition to 
scattering, we have to consider a broad range of wavelengths. 
Therefore, we have to average the mean free path over all fre- 
quencies using a proper weighting function. The proper weight- 
ing function is the energy probability distribution leaving an 
interaction (absorption/reemission or scattering) event which is 
proportional to the emission coefficient r\ v . This emission coeffi- 
cient obeys the following integral equation 
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(13) 



is the wavelength dependent mass scattering coeffi- 



13 1 can be solved iteratively. When an guess of rj v 



where cr v 
cient. Eq. 

is available, an improved estimate can be obtained by substitut- 
ing this guess into the right hand side of Eq. ( 13 1. By repeating 
this procedure several times the solution is easily obtained. A 
good initial guess in most cases is to ignore scattering and take 

T]y = KydBy/dT. 

Usually the procedure described above converges in a few 
iterations. When r\ v is computed (d) and (d 2 ^j are given by 
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Note that when there is no scattering, i.e. Xv — K v, Eqs. ( 12J15 1 
reduce to the Rosseland diffusion coefficient given by Eq 



L2T] 
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The formulae above hold for the case of isotropic scattering, 
which we consider in this paper. However, simple corrections 
can be made for anisotropic scattering. As shown by Ishimaru 
( |1978| l in the diffusion domain similar results are obtained when 
cr v and the asymmetry parameter g are both replaced by cr' v and 
g' given by 

<F (16) 



Table 1. The parameters for the model setup (as used in Eq. 17 1 



cr'„ 



1 - 



The simplest form is by taking the solution for isotropic scatter- 
ing, i.e. g' = 0, with cf v = {\- g)cr v . 

2.3. Description of the code: MCMax 

The radiative transfer code MCMax is based on the Monte Carlo 
method outlined by Bjorkman & Wood]( |2"001| ) to solve the tem- 
perature structure in the disk. The observables are obtained by 
analytic integration of the computed source function in order to 
avoid noise in the observables. The partial diffusion approxima- 
tion as well as the modified random walk as discussed above are 
implemented as options. The code solves the radiative transfer 
in 3 dimensions but it is assumed that the geometry is axisym- 
metric. Sophisticated regridding algorithms make sure that the 
optical depth to the stellar and local radiation field is properly 
sampled so that the temperature structures of extreme density 
distributions can be solved accurately. 

The scattering of radiation by the dust grains is treated in 
a complete sense. This means that the full angle-dependent 
Mueller matrix is used for each scattering event. In this way, 
also polarization maps of the disks can be obtained. A compari- 
son of the results for anisotropic scattering with other codes will 
be presented in Pinte et al. ( submitted] !. In this study we assume 
isotropic scattering. 

In addition to solving the radiative transfer equations subject 
to the constraint of radiative equilibrium, the MCMax code can 
also be used to solve the vertical scale height of the disk self 
consistently by iterating the radiative transfer computations and 
subsequently solving vertical hydrostatic equilibrium equations 
as outlined in section [ 



2.4. Disk setup 

For the density structure of the disk we follow partly the bench- 
mark setup as defined by |Pascucci et al. (2004). However, we 
modify the setup to more massive disks with higher optical 
depths. We deviate from the |Pascucci et al. (2004) setup on three 
important points. 



1. 



The surface density of the disk, S, as chosen by |Pascucci| 
|et al.| ( |2004|) increases with the distance to the central star 
fl as E oc~fl° 125 . This implies that the vertical optical depth 
through the disk increases with distance, putting the most 
optically thick regions far away from the star. We consider 
a radial dependence of the surface density which is a power 
law with £ oc R l (which is closer to the theoretically pre- 
dicted value, see |Dullemond et aL 2007). 
The outer disk radius in the Pascucci et al. ( 2004[ ) setup is 
chosen very large, fl ou t = 1000 AU, while protoplanetary 
disks are thought to have smaller sizes. We therefore adopt a 
value of fl ou t = 200 AU. 

For the vertical scale height of the disk we choose a param- 
eterization which is close to the case of vertical hydrostatic 
equilibrium. In section 3.4 we will show the results for the 



Symbol 


Meaning 


Value 


M* 


Stellar mass 


2.5 M Q 


R* 


Stellar radius 


2 Rq 


T* 


Stellar effective temperature 


10000 K 


Rm 


Inner disk radius 


0.75 AU 


R-uut 


Outer disk radius 


200 AU 


h 


Scale height at 1 AU 


0.03 AU 


P 


Power law for the scale height 


1.3 


a 


Grain radius 


0.12/mi 


Pg 


Grain density 


3.6g/cm 3 


M dust 


Dust mass in the disk 


10- 6 , 10- 3 M o 



For the scale height of the disk we also take a powerlaw with 
distance, h(R) oc RP, where p is determined to be roughly con- 
sistent with that of a disk in vertical hydrostatic equilibrium. The 
density structure we take is for R m < R < R out given by 



M dus ,(lAU)" 



27r -\fn h x (fl out 

m = Miatj/ '■ 



Rm) 



.R-l-P e -(z/h(R)) 2 



(17) 
(18) 



Here p(fl, z) is the dust density in the disk, Md us t is the total dust 
mass, A is the distance to the center of the star (in cylindrical 
coordinates), z is the height above the midplane, p is the power 
of the scaleheight powerlaw, hi is the scaleheight of the disk at 
1 AU, and R m and R oul are the inner and outer radius of the disk. 
For the values of the parameters see Table [T] The radius of the 
inner edge of the disk is chosen to correspond roughly to the dust 
condensation temperature. 
By integrating Eq. 



17 i through the midplane (z = 0) from 



R - Ri n to flout we can compute the radial optical depth 



* radial, 
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Xv 

(19) 

This corresponds, for the parameters of the disks we have cho- 
sen, to an optical depth at 550 nm of T ra di a i,55o nm = 3.5 • 10 6 (3.5 • 
10 3 ) for the model with M dust = 10~ 3 (lO" i5 )M . 

We take a very simplistic description of the grain properties. 
The dust grains in the disk are homogeneous spherical particles 
with a radius of 0.12/zm composed of so-called 'astronomical 
silicate' (Draine & Lee 1984[ ). The scattering by the grains is 
assumed to be isotropic. We would like to note here that al- 
though the grains have in principle a non-isotropic phase func- 
tion, the errors on the temperature structure and emerging spec- 
tra made by assuming isotropic scattering are generally small. 
So although both codes used for the calculations can in principle 
handle anisotropic scattering, we choose to use isotropic scatter- 
ing for the computations in this paper. 

2.5. Vertical structure 

At first we will discuss the results using the above parameterized 
version of the density structure. In section 3.4 we solve the verti- 



iterated solution of the vertical density distribution. 



cal structure using hydrostatic equilibrium computations and it- 
erate on the vertical structure until it sufficiently converges. For 
these computations we kept the surface density fixed at £ oc R l . 
The vertical density profile is then set by the balance between 
pressure and gravity and is given by the solution of the follow- 



M. Min et al.: Radiative transfer in very optically thick circumstellar disks 



5 



ing differential equation (see e.g. |Dullemond| |2002| l. 



dP(R,z) ,„ GM* 

= -p(R,Z) —r- z, 



dz 



(20) 



where G is the gravitational constant, and the pressure P = 
kspT /(jj.m u ) with kg the Boltzmann constant, the mean molec- 
ular weight fi = 2.3 (for a H2, He mixture) and m u is the pro- 
ton mass. The procedure is now to first compute the temper- 
ature structure for an initial guess of the density distribution. 
These temperatures and their derivatives can then be used to 



solve Eq. (20 1, and this can be iterated several times until the 
convergence criterion is met. 

In order to judge if the vertical structure is sufficiently con- 
verged, a convergence criterion is needed. For this we used the 
temperature structure of the disk. The error in the temperature 
structure can be computed using the statistics from the Monte- 
Carlo procedure. This allows us to compute the difference be- 
tween two iterations as compared to the statistical errors. When 
the average difference between the temperature structures of two 
consecutive iterations is less than 3cr we conclude that the model 
is converged. 




Figure 2. Same as the right panel of Fig.[T|but now comparing 
a model using 10 6 photon packages with the reference model 
using 10 8 photon packages. 



3. Results 

3.1. Results with partial diffusion approximation 

In this section we will compare the results with and without the 



partial diffusion approximation (PDA) as outlined in Section 2.2 
We will compare the results of the runs with 10 5 photon pack- 
ages with and without radiative diffusion with those obtained 
using 10 8 photon packages without the partial diffusion approxi- 
mation or modified random walk procedure. We will refer to the 
latter model as the reference model. We have used the PDA ev- 
erywhere where in the 10 s photon model the number of photons 
visiting a cell is below 30 or the relative error in the temperature 
as derived from the photon statistics is higher than 5%. The rela- 
tive difference in the computed temperature structure compared 
to the run with 10 8 photon packages is shown in Fig. [I] In this 
figure negative errors correspond to the case where the temper- 
ature using the PDA is lower than that of the reference model, 
while positive errors represent higher temperatures. It is clear 
that the PDA improves the agreement with the high photon run 
considerably, although there are still regions with relatively large 
errors. These errors are partly caused by errors in the tempera- 
ture at the boundary of the region where the PDA is used. Most 
of the diffusion occurs vertically. Since the cells at the bound- 
ary are used to set the boundary conditions of the PDA, an error 
in the temperature in these cells propagates all the way to the 
midplane, causing the vertical lines of large errors seen in Fig. T] 

The largest errors in the derived temperatures without the 
use of the PDA are 100% (which lies outside the limits of errors 
shown in Fig. fTJ. This is in the regions of the disk where there 
are no photon packages in the 10 5 photon run. By using the PDA 
we can reduce this error to 18%. 

The errors when using only 10 s photon packages and the 
PDA are still considerable. This is because by using such low 
number of photons, the region where the diffusion approxima- 
tion is employed is large. When we increase the number of pho- 
ton packages to 10 6 we can decrease the size of this region sig- 
nificantly (see Fig.|2|. The maximum error in this case is ~ 12%, 
but only in a very small region of the disk. This model still runs 
in a few minutes on a normal desktop computer and because of 
the high gain in accuracy seems to be the best tradeoff between 
speed and accuracy for most cases. 



The diffusion approximation is only valid in regions with 
very high optical depths. The above errors are caused by the 
Monte Carlo noise and the intrinsic errors caused by the as- 
sumption entering the PDA. To study the size of the intrinsic 
errors we have used the partial diffusion approximation to com- 
pute the temperatures in the 10 8 photon run in regions that have 
less than 30000 or 3000 photons. These regions are comparable 
to the regions where the PDA was used for the 10 5 and the 10 6 
photon runs respectively. However, now the temperatures at the 
boundary of the diffusion region are determined to a very high 
accuracy. Therefore, we can now see the intrinsic error made by 
the PDA. For the region corresponding to the 10 5 photon run 
(limit of 30000 photons) the largest relative error in the temper- 
ature compared to the reference model is still ~ 18% (see Fig. [3]>. 
For the region corresponding to the 10 6 photon run (limit 3000 
photons) the largest relative error is ~ 8.5%. 



As mentioned in section 2.2. 1 the PDA always tends to over- 
estimate the temperature. In principle this causes violation of 
energy conservation. When the region where the PDA is used 
is sufficiently shielded from the observer, this is no problem. 
However, one can think of exotic density distributions where 
large, moderately optically thick regions are shielded from stel- 
lar radiation but not from the observer. In these cases, the er- 
rors in the temperature distribution can be high and will result 
in errors in the computed spectra and images when these are 
obtained from the obtained temperature distribution. This will 
be visible because in these cases the integrated total luminosity 
will be higher than the stellar luminosity. Since the PDA is only 
used after the Monte Carlo run, observables directly obtained 
from escaping photons will not be influenced. For ray-tracing 
the observables, like spectra, images and visibilities, the noise 
in the temperature structure usually provides no problems. Thus, 
when ray-tracing observables it is better to use the temperature 
as obtained before using the PDA. For determining the verti- 
cal scale height it is required that the temperature distribution 
is smooth and that the temperature is well defined throughout 
the entire disk. Therefore, temperatures obtained by the PDA 
are used when computing the vertical scaleheight of the disk (or 
other dynamical or geometrical properties of the disk that de- 
pend on temperature and require a smooth temperature profile). 
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Figure 1. The relative error in the temperature computed by the model using 10 5 photon packages as compared to the reference 
model using 10 8 photon packages as a function of position in the disk. The left panel shows the model without the radiative 
diffusion approximation, the right panel with the radiative diffusion approximation. The black contour line encloses the region 
where the number of photon packages entering a cell is below 30, which corresponds to an error in the derived temperature of ~ 5%. 
The dust mass in these models is 10~ 3 Mq. 





Figure 3. The relative error in the temperature computed by the model using 10 8 photon packages with PDA for the region with 
'low' photon counts as compared to the reference model using 10 8 photon packages as a function of position in the disk. Everywhere 
in the disk where less than 30000 photon packages (left panel) or 3000 photon packages (right panel) entered we used the PDA. 
This region is enclosed by the black contour line. The dust mass in these models is 10~ 3 Mq. 



3.2. Results with modified random walk approximation 

In this section we study the accuracy of the Modified Random 
Walk (MRW) procedure. We have run the radiative transfer code 
using 10 8 photon packages with and without the MRW. We used 
y — 1 and y = 10. The difference in execution time is shown 
in Table [2] In this table we present the average runtime per pho- 
ton package for the various models with and without the modi- 
fied random walk (MRW) method relative to the 10~ 6 Mq model 
without the MRW (which was 0.22 ms in our case on a nor- 
mal desktop computer). For this comparison we have also added 
disks with extreme mass (M^ ust = 10" 2 Mq). Although such 
disks are too massive to be stable in reality, it represents the re- 
sults one would get for more extreme density distributions (for 
example a disk with the same mass but with a steeper power law 
index of the surface density). It is clear that the MRW speeds up 
the model for massive disks. 



We compare the two cases with different values for y again 
with the reference model using 10 8 photon packages. The results 
are shown in Fig. |4] In the left panel we show the case for y = 1 . 
The relative error in the temperature in this case varies between 
-17% and +2%. The errors are concentrated in the innermost 
regions. The second case, with y = 10, is shown in the right 
panel of Fig.[4]and has comparable errors. In this case the errors 
vary between -15% and +10%. More importantly, in this case 
the region where the errors occur is much smaller. 

We conclude from the above that the temperatures are, on 
average, systematically underestimated by the MRW. Part of the 
reason for this is that for a small circumscribing sphere, the dif- 
fusion approximation used to derive the random walk procedure 
tends to underestimate the time it takes to leave this sphere. 
When the radius of the sphere is increased the approximation 
gets better. We can do this by increasing y which ensures that 
the MRW is only activated when the circumscribing sphere is 
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Figure 4. The relative error in the temperature using 10 8 photon packages and the MRW compared to the reference model using 10 8 
photon packages as a function of position in the disk. The left panel shows the case where y = 1, the right panel shows the case for 
y = 10. The dust mass in these models is 10~ 3 Mq. 



Table 2. The relative average runtime per photon package for the 
various models. 



M dusl 
[Mq] 


Relative runtime per photon package 
without MRW with MRW 
7=1 y = 10 


io- b 


1 


0.73 


0.73 


io- 3 


10.2 


1.27 


2.05 


io- 2 


58.6 


1.82 


5.23 



large enough. For y — 10 we find that the errors are below 15% 
everywhere in the disk, and that the region with the largest er- 
rors is very small and located near the midplane of the disk. A 
second reason for the errors on the temperature for low values 
of y is that the direction of propagation of the photon package 
after leaving the random walk sphere is not well defined. When 
the radius of the sphere is increased, the travel direction at the 
edge becomes less important and thus the temperature estimate 
becomes more accurate. 

There is no effect of the errors on the temperature structure 
on the observables. This is because the regions where the MRW 
is activated are sufficiently shielded from the observer. However, 
note that in contrast to the PDA, the MRW is used during the 
Monte Carlo procedure and thus, in principle, also influences the 
observables directly obtained from escaping Monte Carlo pho- 
ton packages. These errors will be discussed below. 

The MRW makes the dependence of the computation time 
on the mass of the disk much weaker. In addition to the cases 
presented in this paper, we have performed computations for ex- 
treme density distributions (e.g. with X oc r~ 2 ) where the de- 
crease in computation time due to the random walk procedure is 
up to a few orders of magnitude. 

3.3. Observables: the effect on the SEDs 

To test the effects of the MRW on the predicted observables we 
have considered the spectral energy distribution (SED), visibility 
curves and images. We will discuss here only the results on the 
SEDs. For the case where we solve the hydrostatic equilibrium 
(section |3.4|i we will also show the resulting visibility curves. 



Note that the temperature computed using the PDA is not used in 
the ray-tracing procedure for computing the observables. Instead 
we use the temperature directly derived from the Monte Carlo 
run including the MRW procedure. 

For the SEDs we looked at the entire spectrum from 0.1 jum 
up to 1 cm. The resulting errors for the high mass disk seen under 
an inclination of 45° are plotted in Fig. [3] For the low mass disk 
the effects of using the random walk procedure is negligible, as it 
is only used in a very small region of the disk. For the high mass 
disk the effects of using the MRW with y = 1 can be as large as 
1% for a disk seen under moderate inclinations. For inclinations 
close to edge on (i.e. inclination higher than 80°), the differences 
can be a few (up to 4) percent (not shown). For y = 10 the error 
made using the MRW is much smaller. For moderate inclinations 
it is always lower than 0.2%. For nearly edge on disks the error 
can be up to 0.5% (not shown). The increasing error at the long 
wavelength side of the SED is caused by the fact that at these 
wavelengths one starts to see through the disk, so the temper- 
atures in the midplane are more important. The reason that the 
errors are much smaller in the case of y — 10 is that, although 
the errors on the temperature are of the same order, the region 
where this error is made is much smaller and deeper hidden in 
the disk. We conclude that the effect of the MRW procedure on 
the SED is very small. 

3.4. Solving the vertical density distribution 

The vertical structure of the disk is mostly set by the midplane 
temperature, and therefore it is important to study the effects of 
the PDA and MRW approximations on the iterated vertical den- 
sity distribution. We expect the effects of the errors of the tem- 
perature to translate into somewhat smaller errors in the vertical 
structure since the scale height of the disk is roughly propor- 
tional to y/T. In this section we first analyse how the obtained 
structure of the disk is influenced, and after this evaluate what 
the consequences are for simulated observations. As a conver- 
gence criterion we use the difference between two iterations in 
terms of standard deviations, cr, integrated over the entire disk 
where we use the statistical error in the temperatures to compute 
a normalized standard deviation. 
We consider here three cases: 
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Figure 5. The relative differences between the SEDs computed 
with the random walk procedure using two different values of 
y as compared to the reference model without the random walk 
procedure. The results shown are for the high mass disk under 
an inclination of 45°. 



Case 1: Using 10 5 photon packages, applying the modified ran- 
dom walk (y = 10), and partial diffusion approximation (us- 
ing a minimum number of photons of 30). 

Case 2 Using 10 6 photon packages, applying the modified ran- 
dom walk (y = 10), and partial diffusion approximation (us- 
ing a minimum number of photons of 30). 

Case 3 Using 2 • 10 7 photon packages, without applying the 
modified random walk. We do apply the partial diffusion 
approximation but with a minimum number of photons of 
only 1 (the partial diffusion approximation is only used to as- 
sure that no cells exist with a zero temperature, which would 
cause the vertical scale height to collapse). 

Thus, here Case 3 is our reference model. 

For the low mass disk (Md us t = 10~ 6 Mq) the results of the 
vertical height of the disk is shown in Fig. [6] In this figure we 
show the height of the surface of the disk where an optical depth 
of unity is reached in the direction perpendicular to the midplane 
at a wavelength of 0.55 yum. In the first case (using 10 5 photon 
packages) the diffusion region used for the PDA is quite large, 
causing an overestimate of the midplane temperature and there- 
fore an overestimate of the height of the disk above the diffu- 
sion region. This region lies just behind the inner rim. By using 
10 6 photon packages (Case 2) the diffusion region is sufficiently 
small to overcome this problem, and no significant differences 
are observed between Case 2 and 3. In all cases convergence 
was reached within 9 iterations. 

For the higher mass disk it is much harder to reach conver- 
gence. For this disk mass, for both Case 2 and 3 periodic fluctu- 
ations are observed with a period of 7 to 8 iterations (see Fig. [7}. 
These fluctuations are caused by instabilities in the disk struc- 
ture as described by|D'Alessio et al.| ( [l^9^ ; puiTemond| ( |2000l i; 
Watanabe & Lin ( 2008} and cause the convergence criterion to 
fail. As can be seen from Fig. [7] the solution is periodic, imply- 
ing that the procedure is converged but to an oscillatory solution. 
The question currently remains if in reality these fluctuations are 
suppressed by viscous heating, turbulence or hydrodynamical ef- 
fects causing the disk to respond slower than the heating and 
cooling timescale. For Case 1 (10 5 photon packages) these fluc- 
tuations are not observed and convergence is reached in 7 itera- 
tions (see Pig. [HI. However, the resulting vertical structure is too 



Figure 6. The vertical surface height of the low mass, i.e. Md ust = 
1O~ 6 M0, disk. The ordinate is the height of the surface of the 
disk where an optical depth of unity is reached when looking 
from the top at a wavelength of 0.55 fim divided by the radius, 
R. All three cases described in the text are displayed. The solid 
black line represents Case 1, the dashed line Case 2, and the solid 
grey line Case 3. 




10 15 
iteration number 



Figure 7. The differences in terms of cr (see text) between the 
temperature structures computed in the 26th iteration and all pre- 
vious iterations. Clearly, the structure as in the 26th iteration is 
close to that already computed in the 19th and 12th iteration. 
This oscillatory behavior is caused by 'waves' propagating along 
the surface of the disk. 



high in the region just behind the inner rim. This is, again, caused 
by the overestimate of the temperature at the midplane due to the 
PDA. We conclude that for both cases using 10 6 photon pack- 
ages is sufficient to obtain a reliable vertical density profile. In 
this case the diffusion region, where the PDA is employed, is 
deep enough in the disk to ensure the diffusion approximation to 
be valid and not cause large errors. 

We have computed the spectra and visibility curves of the 
resulting density and temperature structures. The spectra of the 
iterated disks (case 1) together with those of the corresponding 
parametrized disks are shown for disks viewed at an inclination 
of 45° in Fig. [9] In order to judge the effects of the overestimate 
of the disk height just behind the inner rim we show the differ- 
ences in the SEDs between the different cases we considered in 
Fig. 



10 for a disk viewed at an inclination of 15° and 85°. For 



the low mass case it is clear that the largest error is made in the 
l-10yum region. The emission in this spectral region predomi- 



M. Min et al.: Radiative transfer in very optically thick circumstellar disks 



9 




Figure 8. Same as Fig. [fijbut for a dust mass of 10~ 3 Mq. The 
solid black line represents Case 1, the dashed line Case 2, and 
the solid grey line Case 3. The curves for Case 2 and 3 represent 
oscillatory solutions (see main text). 



nantly comes from those radii in the disk where the height of 
the disk is increased because of the use of the PDA (see Fig. [6|, 
causing these regions to emit more efficiently when the PDA is 
employed. The errors at long wavelengths, as seen in the high 
mass case, are caused by the use of the MRW, as discussed above 
for the non-iterated case. In all cases we studied the differences 
are smaller than 5%, independent of disk inclination. Therefore, 
we conclude that for modeling SED observations, it suffices to 
use the very fast, i.e. 10 s photon package, method with the PDA. 
Also, the effects of the oscillations in the vertical structure of the 
massive disk have only a very small impact (less than 5%) on the 
SEDs. 



For the visibilities things can be different. In Fig. 1 1 we show 



the visibility curves in the N-band, as would be observed by e.g. 
the Mid-Infrared Interferometric Instrument (MIDI) at the Very 
Large Telescope Interferometer (VLTI), for the low mass disk at 
a baseline of 100 meters (the disk is set at a distance of 150pc). 
It is already clear that the overestimate in the vertical height for 
Case 1 results in a slightly different visibility profile. For the high 
mass disk we have no convergence of the vertical structure. For 
Case 1 we get very different results as compared to either Case 
2 or 3 which give visibility curves that are very close together 
(not shown). In Fig. 12 we show the visibility curves for Case 2 
after different numbers of iterations. It is clear that the oscillatory 
behavior seen in the surface height of the disk can be observed 
in the visibility curves. 

The oscillations found in the vertical structure of the disk do 
not appear to be caused by numerical problems. The waves orig- 
inate in the region that just emerges from the shadow of the inner 
rim, and is thus again directly illuminated by the central star. In 
our case this is around 8-9 AU. From there onwards they travel 
in until they fade out inside the shadow of the inner rim. The 
fact that it doesn't seem to be a numerical artifact means that 
if there are no processes to stabilize the vertical structure (like 



those discussed by D'Alessio et al. 1999} Dullemond 2000), 



the oscillations in the visibilities mentioned above could in prin- 
ciple be observed. The period with which these oscillations oc- 
cur can be estimated from the time it takes for a sound wave to 
propagate from the midplane to the surface of the disk. Around 
8 AU, where the waves originate, this timescale is ~ 2 years. 
The period of the oscillations is around 7 iterations and thus in 
time is estimated to be roughly 7 times the vertical timescale, i.e. 




X (fJ.m) 



Figure 10. The relative differences between the SEDs computed 
for the three different cases defined in the text. The upper panel 
shows the results for the dust mass of 10~ 6 Mq, the lower for a 
dust mass of 10~ 3 Mq. 




X [^m] 

Figure 11. The visibility curves for the low mass, i.e. Md ust 
10 Mq, disk and a baseline of 100 meter. The distance to the 
object is 150pc. All three cases described in the text are shown. 
The solid black line represents Case 1, the dashed line Case 2, 
and the solid grey line Case 3. 



14 years. We should note that the timescale for a wave to travel 
from the midplane to the surface of the disk is a function of the 
distance to the central star. For example, at 3 AU the timescale 
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Figure 12. The visibility curves for the high mass, i.e. Md ust = 
1(T 3 Mq, disk and a baseline of 100 meter. The distance to the 
object is 150 pc. The visibilities computed after three different 
numbers of iterations for Case 2 are shown. The solid black line 
is after 19 iterations, the dashed line after 23 iterations, and the 
solid grey line after 26 iterations. 

is only ~ 5 months. It is hard to estimate from simple consid- 
erations if this will have an increasing or a decreasing effect on 
the strength of the waves. Computations using radiation hydro- 
dynamics have to be performed to study these effects in more 
detail. 

4. Recommendations for use and implementation 

In this section we summarize the results discussed above and 
make recommendations for use and implementation of the PDA 
and the MRW procedures for various cases. 

4. 1 . Passively heated disks 

For radiative transfer through a fixed density structure the rec- 
ommended setup for computations depends on the observables 
one wants to obtain. For computing SEDs or optical/infrared im- 



ages it suffices to use only the MRW procedure, the PDA is then 
not required. Since the errors made by the MRW procedure are 
localized strongly in the high density regions, one can use a low 
value of y (e.g. ~ 1-5). Errors are then smaller than 1 .5%. When 
also millimeter images or visibilities are required a higher value 
of y is recommended (e.g. ~ 10). 

When accurate temperatures are required, e.g. for computa- 
tions of chemical reactions or of the vertical structure, the PDA 
can be switched on. The number of photon packages used and 
the limit when the PDA is used can be tuned to the desired ac- 
curacy and computational speed. One needs to keep in mind that 
the diffusion approximation itself is not without error, so the use 
of a sufficient number of photon packages is advised. In order to 
keep the statistical error in the temperature structure below ~ 5% 
a lower limit of the number of photon packages per cell has to 
be set to approximately 30. 

There are roughly two types of situations which require dif- 
ferent considerations for the limits of the diffusion approxima- 
tion and the number of photons needed which we will refer to as 
local accuracy and global accuracy. 

4.1.1. Local accuracy: 

The first, and most difficult case is when the temperature is 
needed at each location in the disk with high accuracy, e.g. 5%. 
In this case, the PDA has to be restricted to a very small region 
deep inside the disk to decrease the intrinsic error of the diffusion 
approximation. In order to do so, a sufficient number of photon 
packages (e.g. ~ 10 7 ) have to be used, and the limiting number 
of packages has to be set to ~ 30 to avoid statistical errors. 

4.1.2. Global accuracy: 

The second situation is when it is most important that the tem- 
perature structure is not too noisy, i.e. the fluctuations from cell 
to cell need to be small. This can be the case when one wants 
to look at more global disk characteristics in which the tempera- 
ture plays a role like, for example, when solving for the vertical 
structure of the disk. In this case, a larger error in the tempera- 
ture can be allowed in favor of a more stable and smooth tem- 
perature structure. Consequently, fewer photon packages can be 
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used (e.g. ~ 10 6 ) and the limit for using the PDA can be set to 
~ 30. Although the errors on the temperature structure can in 
this case be about 10%, the noise level is much lower due to the 
smoothing caused by the PDA. Also, in this case the y for the 
MRW procedure can be set to a relatively low value (e.g. y — 1) 
since the temperatures in the region where the MRW is actually 
employed will be largely overwritten by the PDA. 

4.2. Viscously heated disks 

Although in this paper we have not considered the case of vis- 
cous heating, the MRW procedure does allow for complete trans- 
port of the energy released by viscous heating through the disk. 
When the MRW is not employed it is not feasible to transport 
a significant number of photon packages from the midplane, 
where the energy is released, to the surface layers. The compu- 
tations we performed on the high mass disk showed that without 
the MRW procedure the computation time for a photon package 
to reach the surface layers was roughly a minute on a normal 
desktop computer. By using the MRW procedure with y = 1 this 
was reduced, in our case, to less than a second. This still does not 
allow for a very large number of photon packages, but at least 
several thousands can be used. The energy released by viscous 
heating ensures that the regions shielded from stellar radiation 
also receive sufficient photon packages to determine an accurate 
temperature. Thus in this case the PDA can be turned off. More 
details on this subject are beyond the scope of this paper and will 
be left for a future study. 



5. Conclusions 

We have presented two efficient implementations of the diffusion 
approximation in the high density regions of massive protoplan- 
etary disks: the partial diffusion approximation and the modified 
random walk procedure. 

First, the partial diffusion approximation (PDA) can be 
used to increase the accuracy of the temperature structure in 
the highly obscured regions. The PDA is mainly useful when 
the vertical scaleheight of the disk has to be computed self- 
consistently. In order to converge the vertical scaleheight com- 
putations it is important to have a temperature structure that is 
stable throughout the disk. The partial diffusion approximation 
is therefore ideal to reduce the photon noise, inherent to Monte- 
Carlo radiative transfer computations, in the midplane regions. 

Second, the modified random walk (MRW) procedure in- 
creases the computational speed significantly, while errors are 
small. When the modified random walk procedure is not em- 
ployed the computation time required is a strong function of the 
mass of the disk. When it is employed, this dependence is much 
weaker allowing for computations of disks with much more ex- 
treme optical depths, like for example in disks with steeper den- 
sity gradients. 

The combination of the above two methods allows one to 
perform computations for massive protoplanetary disks with 
high accuracy and within reasonable time. The modified random 
walk procedure ensures that enough photons can be emitted in 
order to get reasonable photon counts in most of the disk. The 
partial diffusion approximation can then be employed to correct 
the temperatures in those regions where errors are still large due 
to poor statistics. 

We have performed computations to iteratively compute the 
vertical scaleheight of the disk. For this we have studied care- 
fully the effects of both approximations mentioned above. The 



effects of using the MRW (with y = 10) on the vertical structure 
are very small, and have no impact on the observables. When 
very few photon packages are used, the region where one has to 
apply the PDA is large. This causes overestimates of the temper- 
atures in a relatively large region of the disk, resulting in a local 
overestimate of the vertical height. Although this has no effect 
on the SED, the visibilities obtained with interferometric mea- 
surements are affected. We conclude that for computations of 
spatially unresolved observations, only a small number of pho- 
ton packages is enough. When spatially resolved observations 
are to be simulated, one needs to resort to higher numbers of 
photon packages in order to shrink the region where the PDA 
has to be applied. 

We also find that the vertical structure of massive disks con- 
verges to an oscillatory solution. 'Waves' traveling the surface 
of the disk cause the density structure to fluctuate. These waves 
originate at the location where the disk emerges from the shadow 
of the inner rim and from there travel inwards. When the damp- 
ing mechanisms for these waves are not efficient, this can be 
important when considering spatially resolved computations or 
observations of flaring protoplanetary disks. 
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Appendix A: Benchmark test of MCMax 

The new code MCMax has been tested for accuracy. We did this 
by comparing the computed spectra and temperature structures 
to those obtained by the independently developed radiative trans- 
fer code RADMC (first used in Dulle mond & Dominik"l |2004). 
Both codes have been checked to reproduce the benchmark re- 
sults of P ascucci et al. (2004). In addition to this both codes are 
involved in benchmark computations for higher mass disks by 
Pinte et al. (in prep) and no discrepancies with other codes have 
been found so far. 

We tested the codes MCMax and R ADMC first for the 
benchmark setup of |Pascucci et al.| ( |2004| l. All results where re- 
produced without any problems. The differences with the bench- 
mark SEDs are comparable with those presented by Pascucci 
|et al.| ( |2004[ ) between the different codes used in the bench- 
mark computations. For the temperature structures we compared 
the results with those computed using the RADICAL code and 
found that the differences are in all cases less than 2%. 

As an additional test case we performed computations for 
the 10~ 3 and 10~ 6 solar ma ss dust disks given the density setup 
as described in section 



2.4 



using 10 8 photon packages for both 
codes. By using this number of photon packages we ensure that 
the photon count is very high everywhere in the disk. For this run 
we did not use the random walk or the diffusion approximation 
to ensure that the results are not influenced by the assumptions 
made in these approximations. 

The comparison of the computed temperature structures is 



shown in Fig. A. 1 As the error we plot 



RADMC 



-7l 



MCMax 



TrADMC + ?MCMax 



(A.l) 



For the model with 10~ 6 Mq of dust the differences are below 
6% everywhere in the disk. All differences seem to be caused 
by interpolating between the different spatial grids used in both 
codes. For the model with 10~ 3 Mq of dust, the differences are 
slightly larger. In this case we can distinguish three types of dif- 
ferences: 
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Figure A.l. The relative error in the temperature computed by MCMax and RADMC. The left panel shows the model with a dust 
mass of 10~ 6 Mq, while the right panel shows the model with 1CT 3 Mq of dust. 
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Figure A.2. The resulting SEDs as computed for the two disk masses. The SEDs are shown for two disk inclination angles. The 
solid line gives the resulting spectrum from MCMax while the dotted line shows the result as computed by RADMC. The distance 
to the object is 150 pc. 



1. Those caused by the interpolation of the spatial grids. 

2. Those caused by the fact that the inner boundary in the 
RADMC code is smoothed, while in the MCMax code it is 
sharp. 

3. Errors caused by photon noise and numerical implementa- 
tion of the methods. 

At places in the disk where the temperature gradient is large the 
errors due to interpolation are most important. In these regions, 
these errors can get up to ~ 6%. Differences 2) only play a role 
at the innermost edge of the disk. Here we find errors, which 
are most likely caused by this effect, up to ~ 30%. These errors 
occur only in an extremely small region of the disk, i.e. the in- 
ner 2 • 10~ 4 AU. Note that the temperature of the emitting inner 
rim is the same for both codes, as is the temperature gradient 
expressed in optical depths (which is important for determining 
the inner rim emission). The smoothing of the inner rim as used 
in RADMC, causes a slightly different temperature gradient as 
a function of radius since the optical depth gradient is different. 
The differences occur only inside 5 ■ 10~ 5 AU from the inner rim 



of the dust disk. Errors that cannot be attributed to effects 1 or 2 
are due to numerical implementation of the equations and statis- 
tical noise. These can be considered the true differences caused 
by the differences in the codes. They are typically small, except 
for the midplane region. In this region, the photon count is low 
and the exact implementation of the method becomes very im- 
portant. The true differences are typically ~ 15% close to the 
midplane and smaller than a few percent everywhere else in the 
disk. We thus conclude that overall the two codes agree very well 
on the computed temperature structures. 

The resulting spectral energy distributions for a disk at a 
distance of 150pc at two different inclinations are shown in 
Fig. 



A.2 



The differences between the RADMC and MCMax 
models are small, although for the high inclination cases differ- 
ences can be observed even on a log-log scale. Since the shape 
of the spectra is identical, the most likely cause for this is the 
spatial grid chosen by both codes. 

We conclude that the accuracy of the codes is very good, and 
the results compare very well. 
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